function F = solve_eqm_cf1_inner(param, const, integ, eqm)
%% Counterfactual 1: Baseline

    % fix P_s = 1, e=1.
    P_s     = param.P_s;
    e       = param.e;     

    % read in parameters
    kappa   = param.kappa;
    alpha_m = param.alpha_m; 
    alpha_s = param.alpha_s;
    beta_v  = param.beta_v; 
    beta_m  = param.beta_m;    
    sigma   = param.sigma; 
    f_vH    = param.f_vH;
    f_vF    = param.f_vF;
    mu_E    = param.mu_E;
    sigma_E = param.sigma_E;   
    w       = param.w;
    Q       = param.Q;
    Nf      = param.Nf;
    
    % read in constants    
    gamma   = const.gamma;
    k0      = const.k0;
    C_m     = const.C_m;
    C_v0    = const.C_v0;
    C_x0    = const.C_x0;   
    phi_v   = const.phi_v;
    phi_y   = const.phi_y;
    phi_s   = const.phi_s;
    D_F     = const.D_F;
    lnzQ    = const.lnzQ;
        
    % read in all the integral from previous round
    E_zm0 = integ.E_zm0;
    E_zm1 = integ.E_zm1;
    E_zv0 = integ.E_zv0;
    E_zv1 = integ.E_zv1;
    E_zx0 = integ.E_zx0;
    E_zx1 = integ.E_zx1;
        
    % adjust for level in case integrals too large
    rescale = max(1,10^floor(log10(max(E_zx1))-1));
    E_zx0 = E_zx0/rescale;
    E_zx1 = E_zx1/rescale;
    E_zv0 = E_zv0/(rescale^(1/beta_v));
    E_zv1 = E_zv1/(rescale^(1/beta_v));
    E_zm0 = E_zm0/(rescale^(1/beta_m));        
    E_zm1 = E_zm1/(rescale^(1/beta_m));      

    
%% Solve the equilibrium    

    % Setting
    tol = 1e-12;
    maxIter = 1000;

    % Initial guess 
    D_Hnew = eqm.D_H;
    c_new = eqm.c;    
    
    % Iteration
    err = 1;
    iter = 0;
    while err >= tol && iter < maxIter

        % new guess
        D_H = D_Hnew;
        c = c_new;

        % D(q,E)
        D0 = D_H;
        rv_ads = (D_H/f_vH).^(1/(beta_v-1))./((D_H/f_vH).^(1/(beta_v-1))+(e^sigma*D_F/f_vF).^(1/(beta_v-1)));
        D1 = rv_ads.*D_H + (1-rv_ads)*e^sigma.*D_F;        
        
        % Pi(q)
        Pi0 = D0.^gamma.*(c.^alpha_m*P_s^alpha_s).^((1-sigma)*gamma).*C_x0;
        power = beta_v/(beta_v-1);
        DE = (D_H.^power+(e^sigma*D_F).^power).^(1/power);
        Pi1 = Pi0.*(DE./D_H).^gamma;      
        
        % V(q) and M(q)
        f_v1 = rv_ads.^beta_v*f_vH + (1-rv_ads).^(beta_v)*f_vF;
        C_v1 = (1./(sigma*f_v1.*w)).^(1/beta_v);
        Vbar = C_v0.*Pi0.^(1/beta_v).*E_zv0+rv_ads.*C_v1.*Pi1.^(1/beta_v).*E_zv1;
        V = sum(phi_v.*repmat(Vbar,Q,1),2)';
        M = C_m.*(Pi0.^(1/beta_m).*E_zm0+Pi1.^(1/beta_m).*E_zm1);

        % theta_v(q) and theta_m(q)
        xi = M./V;
        theta_v = 1-exp(-kappa*xi);
        theta_m = (1-exp(-kappa*xi))./xi;
        theta_m(isnan(theta_m))=kappa;       
        theta_m(theta_m==0)=kappa;
        
        % X(q,E)
        X0 = Pi0.*E_zx0;
        X1 = Pi1.*E_zx1;
        X_m = alpha_m*(sigma-1)/sigma*(X0+X1);

        % P(q)^(1-sigma)
        Psig = (X0./D0+rv_ads.*X1./D1);
        %check = Psig.^(1/(1-sigma))./eqm.Psig.^(1/(1-sigma))
        
        % D_m(q)
        TempD = theta_v./M.*c.^(sigma-1).*X_m;
        TempD(isnan(TempD))=0;
        D_m = sum(phi_y.*phi_v.*repmat((TempD)',1,Q),1);

        % X_s
        X_s = eqm.X_s;
        
            % check
            %X = sum(X0 + X1);        
            %B1_implied = 1/sigma*((1-alpha_m-alpha_s)*(sigma-1)+1/beta_v+alpha_m/beta_m+1/gamma)*eqm.X + (alpha_s*(sigma-1)/sigma)*X - X_s;
            rv_rev = rv_ads.*D_H./D1;        
            B1 = sum((1-rv_rev).*X1);
                
        % D_s(q)
        D_s = phi_s.*X_s/sum(phi_s.*Psig);

        % D(q) (new guess)
        D_Hnew = D_m+D_s;

        % c(q) (new guess)
        c_new = (theta_m./V.*sum(phi_y.*phi_v.*repmat(Psig,Q,1),2)').^(1/(1-sigma));           
        
        % update
        D_Hnew = D_H + 0.2*(D_Hnew - D_H);
        c_new = c + 0.2*(c_new - c);
        iter = iter+1;      
        
%             % check 
%             disp('goods market clearing--demand')        
%             [X_s + sum(X_m)]
%             disp('goods market clearing--supply')        
%             [sum(X0)+ sum(rv_rev.*X1)]

        % check convergence
        err1 = max(abs(D_Hnew-D_H)./D_H);
        err2 = max(abs(c_new-c)./c);
        err = err1 + err2;

    end
    
    % Warning
    if iter == maxIter
        msg = 'User Warning: Maximum Iterations Reached.';
        warning(msg)
    end 

    % Adjust back
    Pi0 = Pi0/rescale;
    Pi1 = Pi1/rescale; 
    
    % Save the Export Decision
    ExportCutoffQ = (k0*lnzQ-log(sigma*gamma)+log(kron(ones(Nf,1),Pi1)-kron(ones(Nf,1),Pi0))-mu_E)/sigma_E;
    ExportProbQ = normcdf(ExportCutoffQ);        

    
%% Return

    eqm_cf.Pi0 = Pi0;
    eqm_cf.Pi1 = Pi1;
    eqm_cf.D_H = D_H;
    eqm_cf.D_F = D_F;
    eqm_cf.D0 = D0;
    eqm_cf.D1 = D1;
    eqm_cf.D_s = D_s;
    eqm_cf.D_m = D_m;
    eqm_cf.c = c;
    eqm_cf.Psig = Psig;
    eqm_cf.P = Psig.^(1/(1-sigma));
    eqm_cf.X_m = X_m;
    eqm_cf.X0 = X0;
    eqm_cf.X1 = X1;
    eqm_cf.V = V;
    eqm_cf.Vbar = Vbar;
    eqm_cf.M = M;
    eqm_cf.xi = xi;
    eqm_cf.theta_v = theta_v;
    eqm_cf.theta_m = theta_m;  
    eqm_cf.rv_ads = rv_ads;
    eqm_cf.rv_rev = rv_ads.*D_H./D1; 
    eqm_cf.C_v1 = C_v1;
    eqm_cf.B1 = B1;
    eqm_cf.P_s = P_s;
    eqm_cf.e = e;
    eqm_cf.X_s = X_s;
    eqm_cf.X = sum(X0+X1); 
    eqm_cf.ExportCutoffQ = ExportCutoffQ;
    eqm_cf.ExportProbQ = ExportProbQ;
    eqm_cf.iter_inner = iter;
    eqm_cf.err_inner = err;
    
    F = eqm_cf;

    
end